Logbook Jasper Jonker
Jasper Jonker
2024-16-10
Inleiding
In dit logboek kan je lezen hoe de visualisatie is uitgevoerd op de Transcriptomics van het vak 2.1.2 Genomics & Transcriptomics. In dit vak gingen we bezig met een onderzoek opnieuw doen en kijken welke conclusie wij eruit kunnen halen. Ik ga zelf bezig met het laten zien van de data. Eerst ga ik kijken wat we allemaal kunnen doen en daarna ga ik een vergelijking laten zien.
Hier zie je een aantal library’s die ik heb gebruikt.
library(DESeq2)
library(dplyr)
library(ggplot2)
library(EnhancedVolcano)Mapping en trimmen
mapping is uitgevoerd door ramon om daar meer over te weten lees zijn logbook We hebben 2 soorten cellen, muis en menselijk. De helft van onze samples hebben een muis cellen en de andere helft menselijke cellen. We gaan eerst de mapping doen tegen de muiscellen, die is geindexed (zie janine haar logboek)
cat /students/2024-2025/Thema05/BlaasKanker/Transcriptomics/mouse_cell_SRR.txt | \
parallel 'STAR --runThreadN 6 ' \
'--genomeDir /students/2024-2025/Thema05/BlaasKanker/Transcriptomics/tools/star/index_reference/ ' \
'--readFilesIn /students/2024-2025/Thema05/BlaasKanker/Transcriptomics/output/fastq/{}_1.fastq /students/2024-2025/Thema05/BlaasKanker/Transcriptomics/output/fastq/{}_2.fastq ' \
'--outSAMtype BAM SortedByCoordinate ' \
'--quantMode GeneCounts ' \
'--genomeLoad LoadAndRemove' \
'--limitBAMsortRAM 2000000000 ' \
'--outFileNamePrefix /students/2024-2025/Thema05/BlaasKanker/Transcriptomics/output/mapping/{}_star_'Er is ook op aan geadviseerd om uit eindelijk niet te trimmen
Deseq
Hieronder staat de code voor de verwerking van deseq dit is meer ook een voorbeeld over hoe je het moet gebruiken. Dit is dus van alle groepen van ons onderzoek. Dus tussen NSG (slecht immuunsysteem) muizen en C57BL/6 (goed immuunsysteem) muizen. Elke van deze drie groepen zijn weer onderverdeeld in 3 groepen: Vehicle, Baseline en Entinostat. Ik ga vergelijkingen laten zien van C57BL/6 Baseline tegen C57BL/6 Entinostat. Basline is het Blaaskanker uiteindelijk uit gehaald bij een bepaalde grote en bij Entinostat is er bij die grote Entinostat toegevoegd om te kijken wat het verschil in uiting is
Dit stukje code is van Ramon dus als je er meer uit leg over wilt hebben moet je in zijn logbook kijken.
file.names <- list.files('/students/2024-2025/Thema05/BlaasKanker/Transcriptomics/output/mapping/',
pattern = '*_star_ReadsPerGene.out.tab')
## Function for reading in files
read_sample <- function(file.name) {
## Extract the sample name for naming the column (retaining the 'SRR....' part)
sample.name <- strsplit(file.name, ".", fixed = TRUE)[[1]][1]
sample <- read.table(file.name, header = FALSE, sep="\t",
row.names = NULL, skip = 4)
## Rename the count column
names(sample)[2] <- sample.name
## Return a subset containing the transcript ID and sample name columns
return(sample[c(1, 2)])
}setwd('/students/2024-2025/Thema05/BlaasKanker/Transcriptomics/output/mapping/')
## Read the FIRST sample
counts <- read_sample(file.names[1])
## Read the remaining files and merge the contents
for (file.name in file.names[2:length(file.names)]) {
sample <- read_sample(file.name)
counts <- merge(counts, sample, by = 1)
}
# Set the row names to the transcript IDs
rownames(counts) <- counts$V1
counts <- counts[-1]col_data <- read.csv("/students/2024-2025/Thema05/BlaasKanker/etc/DSEQ_verwerking(Sheet1).csv", sep = ";")
rownames(col_data) <- col_data[,1]
col_data$source_name## [1] "C57BL/6_bladder tumor_Vehicle" "C57BL/6_bladder tumor_Entinostat"
## [3] "C57BL/6_bladder tumor_Vehicle" "C57BL/6_bladder tumor_Vehicle"
## [5] "C57BL/6_bladder tumor_Vehicle" "C57BL/6_bladder tumor_Entinostat"
## [7] "C57BL/6_bladder tumor_Entinostat" "C57BL/6_bladder tumor_Vehicle"
## [9] "C57BL/6_bladder tumor_Baseline" "C57BL/6_bladder tumor_Baseline"
## [11] "C57BL/6_bladder tumor_Entinostat" "C57BL/6_bladder tumor_Baseline"
## [13] "C57BL/6_bladder tumor_Entinostat" "C57BL/6_bladder tumor_Entinostat"
## [15] "C57BL/6_bladder tumor_Vehicle" "C57BL/6_bladder tumor_Baseline"
## [17] "NSG_bladder tumor_Baseline" "NSG_bladder tumor_Baseline"
## [19] "NSG_bladder tumor_Baseline" "NSG_bladder tumor_Baseline"
## [21] "NSG_bladder tumor_Baseline" "NSG_bladder tumor_Entinostat"
## [23] "NSG_bladder tumor_Entinostat" "NSG_bladder tumor_Entinostat"
## [25] "NSG_bladder tumor_Entinostat" "NSG_bladder tumor_Entinostat"
## [27] "NSG_bladder tumor_Vehicle" "NSG_bladder tumor_Vehicle"
## [29] "NSG_bladder tumor_Vehicle" "NSG_bladder tumor_Vehicle"
## [31] "NSG_bladder tumor_Vehicle"Annot count_data
Het annoteren van onze data zorgt ervoor dat de genen de NCBI naamgeving krijgen. Zo is het straks overzichtelijker om informatie te vinden over de genen.
library(dplyr)
library(tidyr)
col_data <- col_data %>%
dplyr::group_by(strain, treatment) %>%
dplyr::mutate(r_num = row_number()) %>%
dplyr::ungroup() %>%
mutate(condition = paste0(strain, "_", treatment, "_r", r_num))
head(col_data)## # A tibble: 6 × 6
## Run source_name strain treatment r_num condition
## <chr> <chr> <chr> <chr> <int> <chr>
## 1 SRR12129014_star_ReadsPerGene C57BL/6_bladde… C57BL… Vehicle 1 C57BL/6_…
## 2 SRR12129015_star_ReadsPerGene C57BL/6_bladde… C57BL… Entinost… 1 C57BL/6_…
## 3 SRR12129016_star_ReadsPerGene C57BL/6_bladde… C57BL… Vehicle 2 C57BL/6_…
## 4 SRR12129017_star_ReadsPerGene C57BL/6_bladde… C57BL… Vehicle 3 C57BL/6_…
## 5 SRR12129018_star_ReadsPerGene C57BL/6_bladde… C57BL… Vehicle 4 C57BL/6_…
## 6 SRR12129019_star_ReadsPerGene C57BL/6_bladde… C57BL… Entinost… 2 C57BL/6_…In de bestanden kon niks gevonden worden over de genen, maar zo is het wel duidelijker waar elke kolom voor staat. Dat moet nu nog toegepast worden in de counts df
for (i in 1:nrow(col_data)) {
idx <- grep(col_data$Run[i], names(counts))
names(counts)[idx] <- col_data$condition[i]
}
head(counts)## C57BL/6_Vehicle_r1 C57BL/6_Entinostat_r1 C57BL/6_Vehicle_r2
## MissingGeneID 5859 4387 7820
## gene-0610005C13Rik 3 0 0
## gene-0610006L08Rik 0 1 0
## gene-0610009E02Rik 0 1 1
## gene-0610009L18Rik 53 29 25
## gene-0610010K14Rik 1024 790 881
## C57BL/6_Vehicle_r3 C57BL/6_Vehicle_r4 C57BL/6_Entinostat_r2
## MissingGeneID 8707 3405 3821
## gene-0610005C13Rik 0 0 0
## gene-0610006L08Rik 0 0 2
## gene-0610009E02Rik 1 0 1
## gene-0610009L18Rik 93 43 34
## gene-0610010K14Rik 1408 897 523
## C57BL/6_Entinostat_r3 C57BL/6_Vehicle_r5 C57BL/6_Baseline_r1
## MissingGeneID 3777 6823 4673
## gene-0610005C13Rik 1 0 0
## gene-0610006L08Rik 0 0 0
## gene-0610009E02Rik 1 4 0
## gene-0610009L18Rik 44 93 25
## gene-0610010K14Rik 492 1304 975
## C57BL/6_Baseline_r2 C57BL/6_Entinostat_r4
## MissingGeneID 7162 2160
## gene-0610005C13Rik 0 1
## gene-0610006L08Rik 0 1
## gene-0610009E02Rik 0 1
## gene-0610009L18Rik 55 15
## gene-0610010K14Rik 1390 419
## C57BL/6_Baseline_r3 C57BL/6_Entinostat_r5
## MissingGeneID 5096 2718
## gene-0610005C13Rik 0 0
## gene-0610006L08Rik 0 0
## gene-0610009E02Rik 0 4
## gene-0610009L18Rik 24 49
## gene-0610010K14Rik 972 676
## C57BL/6_Entinostat_r6 C57BL/6_Vehicle_r6 C57BL/6_Baseline_r4
## MissingGeneID 3870 9562 5163
## gene-0610005C13Rik 0 0 0
## gene-0610006L08Rik 0 0 0
## gene-0610009E02Rik 6 0 2
## gene-0610009L18Rik 34 43 45
## gene-0610010K14Rik 431 964 876
## NSG_Baseline_r1 NSG_Baseline_r2 NSG_Baseline_r3
## MissingGeneID 8453 8045 5803
## gene-0610005C13Rik 4 1 2
## gene-0610006L08Rik 0 0 0
## gene-0610009E02Rik 7 4 7
## gene-0610009L18Rik 38 47 31
## gene-0610010K14Rik 967 1298 664
## NSG_Baseline_r4 NSG_Baseline_r5 NSG_Entinostat_r1
## MissingGeneID 6777 8639 6767
## gene-0610005C13Rik 1 3 1
## gene-0610006L08Rik 1 0 0
## gene-0610009E02Rik 8 10 6
## gene-0610009L18Rik 45 48 95
## gene-0610010K14Rik 1001 1336 1470
## NSG_Entinostat_r2 NSG_Entinostat_r3 NSG_Entinostat_r4
## MissingGeneID 9366 9204 3295
## gene-0610005C13Rik 6 2 0
## gene-0610006L08Rik 0 0 0
## gene-0610009E02Rik 5 6 0
## gene-0610009L18Rik 139 124 64
## gene-0610010K14Rik 1528 1339 387
## NSG_Entinostat_r5 NSG_Vehicle_r1 NSG_Vehicle_r2
## MissingGeneID 12470 5393 7122
## gene-0610005C13Rik 2 2 1
## gene-0610006L08Rik 1 0 0
## gene-0610009E02Rik 13 6 11
## gene-0610009L18Rik 113 41 50
## gene-0610010K14Rik 1393 735 1145
## NSG_Vehicle_r3 NSG_Vehicle_r4 NSG_Vehicle_r5
## MissingGeneID 8495 8854 7730
## gene-0610005C13Rik 2 4 2
## gene-0610006L08Rik 0 0 1
## gene-0610009E02Rik 5 8 4
## gene-0610009L18Rik 122 150 160
## gene-0610010K14Rik 1618 1689 1511Nu zijn de SRR* namen vervangen met waar ze voor staan en is de df duidelijker te lezen, om straks makkelijker de kolommen op te halen die horen bij elk mogelijke variant groepeer ik deze.
print(colnames(counts))## [1] "C57BL/6_Vehicle_r1" "C57BL/6_Entinostat_r1" "C57BL/6_Vehicle_r2"
## [4] "C57BL/6_Vehicle_r3" "C57BL/6_Vehicle_r4" "C57BL/6_Entinostat_r2"
## [7] "C57BL/6_Entinostat_r3" "C57BL/6_Vehicle_r5" "C57BL/6_Baseline_r1"
## [10] "C57BL/6_Baseline_r2" "C57BL/6_Entinostat_r4" "C57BL/6_Baseline_r3"
## [13] "C57BL/6_Entinostat_r5" "C57BL/6_Entinostat_r6" "C57BL/6_Vehicle_r6"
## [16] "C57BL/6_Baseline_r4" "NSG_Baseline_r1" "NSG_Baseline_r2"
## [19] "NSG_Baseline_r3" "NSG_Baseline_r4" "NSG_Baseline_r5"
## [22] "NSG_Entinostat_r1" "NSG_Entinostat_r2" "NSG_Entinostat_r3"
## [25] "NSG_Entinostat_r4" "NSG_Entinostat_r5" "NSG_Vehicle_r1"
## [28] "NSG_Vehicle_r2" "NSG_Vehicle_r3" "NSG_Vehicle_r4"
## [31] "NSG_Vehicle_r5"C57BL_Vehicle <- grep("C57BL/6_Vehicle", names(counts))
C57BL_Entinostat <- grep("C57BL/6_Entinostat", names(counts))
C57BL_Baseline <- grep("C57BL/6_Baseline", names(counts))
nsg_Vehicle <- grep("NSG_Vehicle", names(counts))
nsg_Entinostat <- grep("NSG_Entinostat", names(counts))
nsg_Baseline <- grep("NSG_Baseline", names(counts))Nu kunnen deze variabelen gebruikt worden om de juiste kolommen te selecteren.
Dan kan dit nu samengevoegd worden met DESEQ, wat Janine heeft uitgezocht. Bekijk haar logboek voor meer informatie over de code.
dds <- DESeqDataSetFromMatrix(countData = counts,
colData = col_data,
design = ~ 0 + source_name)
head(dds$source_name)## [1] C57BL/6_bladder tumor_Vehicle C57BL/6_bladder tumor_Entinostat
## [3] C57BL/6_bladder tumor_Vehicle C57BL/6_bladder tumor_Vehicle
## [5] C57BL/6_bladder tumor_Vehicle C57BL/6_bladder tumor_Entinostat
## 6 Levels: C57BL/6_bladder tumor_Baseline ... NSG_bladder tumor_Vehicle# Prefiltering on low gene counts
keep <- rowSums(counts(dds)) >= 10
ddst <- dds[keep,]# Setting factor level
#dds$treatment <- relevel(dds$treatment, ref = "Baseline")
dds$treatment <- factor(dds$treatment, levels = c("Baseline","Entinostat", "Vehicle"), )
dds$treatment## [1] Vehicle Entinostat Vehicle Vehicle Vehicle Entinostat
## [7] Entinostat Vehicle Baseline Baseline Entinostat Baseline
## [13] Entinostat Entinostat Vehicle Baseline Baseline Baseline
## [19] Baseline Baseline Baseline Entinostat Entinostat Entinostat
## [25] Entinostat Entinostat Vehicle Vehicle Vehicle Vehicle
## [31] Vehicle
## Levels: Baseline Entinostat Vehicle# Running deseq
dds <- DESeq(dds)
res <- results(dds)head(res)## log2 fold change (MLE): source_name NSG_bladder tumor_Vehicle vs C57BL/6_bladder tumor_Baseline
## Wald test p-value: source_name NSG_bladder tumor_Vehicle vs C57BL/6_bladder tumor_Baseline
## DataFrame with 6 rows and 6 columns
## baseMean log2FoldChange lfcSE stat pvalue
## <numeric> <numeric> <numeric> <numeric> <numeric>
## MissingGeneID 6213.959413 0.0613458 0.281620 0.2178319 0.82756011
## gene-0610005C13Rik 0.928517 3.1919781 1.353181 2.3588702 0.01833067
## gene-0610006L08Rik 0.325219 0.1850353 3.878567 0.0477071 0.96194965
## gene-0610009E02Rik 3.241600 3.3513189 1.108739 3.0226393 0.00250581
## gene-0610009L18Rik 60.496770 1.1274910 0.398268 2.8309856 0.00464048
## gene-0610010K14Rik 998.035854 -0.0444090 0.294859 -0.1506110 0.88028262
## padj
## <numeric>
## MissingGeneID 0.8980987
## gene-0610005C13Rik 0.0596486
## gene-0610006L08Rik NA
## gene-0610009E02Rik 0.0128493
## gene-0610009L18Rik 0.0208305
## gene-0610010K14Rik 0.9335965Dit is niet wat ik wou hebben of kan gebruiken voor mijn gedeelte. Nu is er een vergelijking gedaan van NSG_bladder tumor_Vehicle tegen C57BL/6_bladder tumor_Baseline. Ik ga een subset maken van alle C57BL/6 baseline en C57BL/6 Entinostat. Deze twee situaties ga ik dus uitwerken.
C57BL_subset <- dds[, c(C57BL_Entinostat,C57BL_Baseline)]
C57BL_subset <- DESeqDataSet(C57BL_subset, design = ~ treatment)
C57BL_subset$treatment <- relevel(C57BL_subset$treatment, ref = "Baseline")
C57BL_subset <- DESeq(C57BL_subset)
resultaat_C57BL <- results(C57BL_subset)head(resultaat_C57BL)## log2 fold change (MLE): treatment Entinostat vs Baseline
## Wald test p-value: treatment Entinostat vs Baseline
## DataFrame with 6 rows and 6 columns
## baseMean log2FoldChange lfcSE stat pvalue
## <numeric> <numeric> <numeric> <numeric> <numeric>
## MissingGeneID 6168.262965 0.17444819 0.261976 0.66589386 0.5054789
## gene-0610005C13Rik 0.405535 1.77823220 3.225919 0.55123283 0.5814741
## gene-0610006L08Rik 0.796704 2.73551203 2.365204 1.15656470 0.2474503
## gene-0610009E02Rik 2.715126 2.96435360 1.375938 2.15442441 0.0312069
## gene-0610009L18Rik 52.611245 0.70508328 0.357148 1.97420351 0.0483586
## gene-0610010K14Rik 1092.089170 0.00181008 0.364426 0.00496695 0.9960370
## padj
## <numeric>
## MissingGeneID 0.6245088
## gene-0610005C13Rik NA
## gene-0610006L08Rik NA
## gene-0610009E02Rik 0.0737484
## gene-0610009L18Rik 0.1045137
## gene-0610010K14Rik 0.9977485Visualization oefen
Ik ga als eerst bezig met alles laten zien en het uitzoeken van mogelijke manieren om het te laten zien. Dit doe ik omdat het niet nodig is om met zijn vieren bezig te gaan met trimmen en deseq. We hebben niet alle tijd, dus moeten we slim bezig gaan met onze tijd.
Pathway Analysis
Hier is een pathway analysis. Hier kan je dus zien als een gen word aangetast en wat voor invloed dit heeft op bepaalde biologische processen.
pathway_pdf is een een pdf met alle informatie
library(pathview)
data(gse16873.d)
pv.out <- pathview(gene.data = gse16873.d[, 1], pathway.id = "04110",
species = "hsa", out.suffix = "gse16873")
Volcano plot oefen
Een volcano plot is handig om te gebruiken, want daarin kan je zien hoe groot een verandering is op de x-as (de lof2FoldChange) en hoe significant een verandering is op de y-as(p-value).
Manier 1
Dit is 1 manier om een volcano plot te maken. Op deze manier wordt het hem denk ik niet. Het word hier nogal moeilijk gemaakt, terwijl het makkelijker kan.
tmp <- readRDS("/students/2024-2025/Thema05/BlaasKanker/Transcriptomics/testData/de_df_for_volcano.rds")
de <- tmp[complete.cases(tmp), ]
de$delabel <- NA
de$diffexpressed <- "NO"
de$delabel[de$diffexpressed != "NO"] <- de$gene_symbol[de$diffexpressed != "NO"]
de$diffexpressed[de$log2FoldChange > 0.6 & de$pvalue < 0.05] <- "UP"
# if log2Foldchange < -0.6 and pvalue < 0.05, set as "DOWN"
de$diffexpressed[de$log2FoldChange < -0.6 & de$pvalue < 0.05] <- "DOWN"
ggplot(data=de, aes(x=log2FoldChange, y=-log10(pvalue), col=diffexpressed, label=delabel)) +
geom_point() +
theme_minimal() +
geom_text() +
geom_vline(xintercept=c(-0.6, 0.6), col="red") +
geom_hline(yintercept=-log10(0.05), col="red") +
scale_color_manual(values=c("blue", "black", "purple"))
Manier 2
Hier zie je ook dat het met Deseq wordt gerund. Hier is het makkelijker om te laten zien. Ook zie je gelijk hoe de deseq werkt en hoe makkelijk je er data uit kan halen.
library(airway)
library(magrittr)
data('airway')
airway$dex %<>% relevel('untrt')
ens <- rownames(airway)library(org.Hs.eg.db)
symbols <- mapIds(org.Hs.eg.db, keys = ens,
column = c('SYMBOL'), keytype = 'ENSEMBL')
symbols <- symbols[!is.na(symbols)]
symbols <- symbols[match(rownames(airway), names(symbols))]
rownames(airway) <- symbols
keep <- !is.na(rownames(airway))
airway <- airway[keep,]library('DESeq2')
dds <- DESeqDataSet(airway, design = ~ cell + dex)
dds <- DESeq(dds, betaPrior=FALSE)
res <- results(dds,
contrast = c('dex','trt','untrt'))
res_1 <- lfcShrink(dds,
contrast = c('dex','trt','untrt'), res=res, type = 'normal') EnhancedVolcano(res_1,
lab = rownames(res_1),
x = 'log2FoldChange',
y = 'pvalue')
Clusterprofiling gsa oefen
Overview of enrichment analysis
library("clusterProfiler")
data(geneList, package="DOSE")
head(geneList)## 4312 8318 10874 55143 55388 991
## 4.572613 4.514594 4.418218 4.144075 3.876258 3.677857gene <- names(geneList)[abs(geneList) > 2]
head(gene)## [1] "4312" "8318" "10874" "55143" "55388" "991"GO classification
ggo <- groupGO(gene = gene,
OrgDb = org.Hs.eg.db,
ont = "CC",
level = 3,
readable = TRUE)
head(ggo)## ID Description Count GeneRatio geneID
## GO:0000133 GO:0000133 polarisome 0 0/207
## GO:0000408 GO:0000408 EKC/KEOPS complex 0 0/207
## GO:0000417 GO:0000417 HIR complex 0 0/207
## GO:0000444 GO:0000444 MIS12/MIND type complex 0 0/207
## GO:0000808 GO:0000808 origin recognition complex 0 0/207
## GO:0000930 GO:0000930 gamma-tubulin complex 0 0/207GO over-representation analysis
Hier kan je zien dat bepaalde genen een overrepresentatie hebben, omdat ze een p.adjust hebben die kleiner is dan 0.05 en de qvalue kleiner is dan 0.05. Dat betekent dat de deze genen meer tot uiting komen dan verwacht wordt. Gene Set Enrichment Analysis
ego <- enrichGO(gene = gene,
universe = names(geneList),
OrgDb = org.Hs.eg.db,
ont = "CC",
pAdjustMethod = "BH",
pvalueCutoff = 0.01,
qvalueCutoff = 0.05,
readable = TRUE)
head(ego)## ID Description GeneRatio
## GO:0005819 GO:0005819 spindle 26/201
## GO:0000775 GO:0000775 chromosome, centromeric region 19/201
## GO:0072686 GO:0072686 mitotic spindle 17/201
## GO:0000779 GO:0000779 condensed chromosome, centromeric region 16/201
## GO:0098687 GO:0098687 chromosomal region 23/201
## GO:0000793 GO:0000793 condensed chromosome 19/201
## BgRatio pvalue p.adjust qvalue
## GO:0005819 332/11884 6.308871e-11 1.886352e-08 1.713356e-08
## GO:0000775 188/11884 3.940180e-10 5.890568e-08 5.350349e-08
## GO:0072686 151/11884 6.423389e-10 6.401978e-08 5.814858e-08
## GO:0000779 138/11884 1.350430e-09 1.009446e-07 9.168708e-08
## GO:0098687 305/11884 1.819474e-09 1.088045e-07 9.882616e-08
## GO:0000793 210/11884 2.580005e-09 1.285703e-07 1.167792e-07
## geneID
## GO:0005819 CDCA8/CDC20/KIF23/CENPE/ASPM/DLGAP5/SKA1/NUSAP1/TPX2/TACC3/NEK2/CDK1/MAD2L1/KIF18A/BIRC5/KIF11/TRAT1/TTK/AURKB/PRC1/KIFC1/KIF18B/KIF20A/AURKA/CCNB1/KIF4A
## GO:0000775 CDCA8/CENPE/NDC80/TOP2A/HJURP/SKA1/NEK2/CENPM/CENPN/ERCC6L/MAD2L1/KIF18A/CDT1/BIRC5/EZH2/TTK/NCAPG/AURKB/CCNB1
## GO:0072686 KIF23/CENPE/ASPM/SKA1/NUSAP1/TPX2/TACC3/CDK1/MAD2L1/KIF18A/KIF11/TRAT1/AURKB/PRC1/KIFC1/KIF18B/AURKA
## GO:0000779 CENPE/NDC80/HJURP/SKA1/NEK2/CENPM/CENPN/ERCC6L/MAD2L1/KIF18A/CDT1/BIRC5/TTK/NCAPG/AURKB/CCNB1
## GO:0098687 CDCA8/CENPE/NDC80/TOP2A/HJURP/SKA1/NEK2/CENPM/RAD51AP1/CENPN/CDK1/ERCC6L/MAD2L1/KIF18A/CDT1/BIRC5/EZH2/TTK/NCAPG/AURKB/CHEK1/CCNB1/MCM5
## GO:0000793 CENPE/NDC80/TOP2A/NCAPH/HJURP/SKA1/NEK2/CENPM/CENPN/ERCC6L/MAD2L1/KIF18A/CDT1/BIRC5/TTK/NCAPG/AURKB/CHEK1/CCNB1
## Count
## GO:0005819 26
## GO:0000775 19
## GO:0072686 17
## GO:0000779 16
## GO:0098687 23
## GO:0000793 19GO Gene Set Enrichment Analysis
Gene Set Enrichment Analysis is een video om Gene Set Enrichment Analysis resultaten uit te leggen hoe werkt.
How to interpret GSEA results and plot is een video die de resultaten uitlegt.
Dus je ziet hier GO:0000775 het meest interessant is door NES met een score van 0.6230073
ego3 <- gseGO(geneList = geneList,
OrgDb = org.Hs.eg.db,
ont = "CC",
minGSSize = 100,
maxGSSize = 500,
pvalueCutoff = 0.05,
verbose = FALSE)
head(ego3)## ID Description setSize
## GO:0000775 GO:0000775 chromosome, centromeric region 188
## GO:0000779 GO:0000779 condensed chromosome, centromeric region 138
## GO:0000776 GO:0000776 kinetochore 130
## GO:0000228 GO:0000228 nuclear chromosome 175
## GO:0098687 GO:0098687 chromosomal region 305
## GO:0000793 GO:0000793 condensed chromosome 210
## enrichmentScore NES pvalue p.adjust qvalue rank
## GO:0000775 0.5970689 2.657162 1e-10 1.62e-09 1.052632e-09 530
## GO:0000779 0.6216009 2.657052 1e-10 1.62e-09 1.052632e-09 798
## GO:0000776 0.6230073 2.649096 1e-10 1.62e-09 1.052632e-09 449
## GO:0000228 0.5875327 2.599594 1e-10 1.62e-09 1.052632e-09 1905
## GO:0098687 0.5419949 2.564564 1e-10 1.62e-09 1.052632e-09 1721
## GO:0000793 0.5507360 2.493279 1e-10 1.62e-09 1.052632e-09 1721
## leading_edge
## GO:0000775 tags=20%, list=4%, signal=19%
## GO:0000779 tags=24%, list=6%, signal=23%
## GO:0000776 tags=21%, list=4%, signal=20%
## GO:0000228 tags=34%, list=15%, signal=29%
## GO:0098687 tags=27%, list=14%, signal=24%
## GO:0000793 tags=28%, list=14%, signal=24%
## core_enrichment
## GO:0000775 55143/1062/10403/7153/55355/220134/4751/79019/55839/54821/4085/81930/81620/332/2146/7272/64151/9212/891/11004/5347/701/11130/79682/57405/10615/79075/2491/11339/3070/9918/1058/699/1063/55055/1051/8970
## GO:0000779 1062/10403/55355/220134/4751/79019/55839/54821/4085/81930/81620/332/7272/64151/9212/891/11004/5347/701/11130/79682/57405/10615/2491/9918/1058/699/1063/55055/1051/79980/9735/23310
## GO:0000776 1062/10403/55355/220134/4751/79019/55839/54821/4085/81930/81620/332/7272/9212/891/11004/5347/701/11130/79682/57405/10615/2491/1058/699/1063/55055
## GO:0000228 8318/55388/7153/23397/4751/9837/332/64151/51659/1111/4174/4171/5347/5888/23594/4998/4175/4173/54962/9918/84296/81611/5111/64785/55506/641/1869/54892/8357/5427/23649/4176/58516/5557/3066/11169/8607/10051/1104/5558/4172/5424/5885/7283/10592/8914/51377/86/5393/142/6839/23212/3014/3619/5425/54107/11177/7273/6119
## GO:0098687 55143/1062/10403/7153/55355/220134/4751/79019/10635/55839/983/54821/4085/81930/81620/332/2146/7272/64151/9212/1111/891/4174/4171/11004/5347/701/11130/79682/57405/10615/79075/2491/5888/4998/11339/3070/4175/4173/2237/9918/1058/699/1063/5111/9401/55055/55506/641/1763/1051/8970/4176/79980/4436/9735/23310/3066/5499/10051/4172/5424/5885/11200/2072/92822/54908/6396/4683/54069/86/3018/1788/142/6839/23028/1786/3014/5905/3619/11335/55166
## GO:0000793 1062/10403/7153/23397/55355/220134/4751/79019/55839/54821/4085/81930/81620/332/7272/64151/9212/1111/891/11004/5347/701/11130/79682/57405/10615/2491/5888/4288/9918/1058/699/1063/55055/641/1051/54892/3148/79980/9735/23310/10051/1104/5885/7283/92822/54908/10592/6396/86/6839/23212/8815/3014/5905/3619/11335/55166p1 <- gseaplot(ego3, geneSetID = 1, by = "runningScore", title = ego3$Description[1])
p2 <- gseaplot(ego3, geneSetID = 1, by = "preranked", title = ego3$Description[1])
cowplot::plot_grid(p1, p2, ncol=1, labels=LETTERS[1:3])
Visualization of functional enrichment result is een handige website met meer mogelijkheiden van laten zien
Er staan een aantal viedeos tussen die zeker handig zijn om te kijken om het beter te bergrijpen.
Visualization C57BL ent tegen Basline
Clusterprofiling gsa
Hier ga ik de resultaten laten zien van de enrichment result. Je kan het op veel mogelijke manieren laten zien, dus laten we eerst maar even de resultaten maken. Hieronder zie je een aantal libraries die ik ga gebruiken bij de visualisatie van de plotjes.
library(clusterProfiler)
library(org.Mm.eg.db)
library(ggplot2)
library(enrichplot)Voor de enrichment result heb je de deseq nodig van je data. Hieronder zie je nog wel even snel de deseq voor treatment Entinostat vs Baseline in de C57BL/6 muizen soort.
head(resultaat_C57BL)## log2 fold change (MLE): treatment Entinostat vs Baseline
## Wald test p-value: treatment Entinostat vs Baseline
## DataFrame with 6 rows and 6 columns
## baseMean log2FoldChange lfcSE stat pvalue
## <numeric> <numeric> <numeric> <numeric> <numeric>
## MissingGeneID 6168.262965 0.17444819 0.261976 0.66589386 0.5054789
## gene-0610005C13Rik 0.405535 1.77823220 3.225919 0.55123283 0.5814741
## gene-0610006L08Rik 0.796704 2.73551203 2.365204 1.15656470 0.2474503
## gene-0610009E02Rik 2.715126 2.96435360 1.375938 2.15442441 0.0312069
## gene-0610009L18Rik 52.611245 0.70508328 0.357148 1.97420351 0.0483586
## gene-0610010K14Rik 1092.089170 0.00181008 0.364426 0.00496695 0.9960370
## padj
## <numeric>
## MissingGeneID 0.6245088
## gene-0610005C13Rik NA
## gene-0610006L08Rik NA
## gene-0610009E02Rik 0.0737484
## gene-0610009L18Rik 0.1045137
## gene-0610010K14Rik 0.9977485Maar eerst moeten wij een genlijst maken want die heeft gseGO nodig om de enrichment russults te doen. De genlijst maken we met behulp van de deseq. Tijdens het maken van de genlijst halen we de NA resultaten van de log2FoldChange eruit. Hieronder zie je de code die er voor gebruikt word. We sorteren het gelijk op volgorde van hoog naar laag.
resultaat_C57BL_2 <- resultaat_C57BL[!is.na(resultaat_C57BL$log2FoldChange),]
gen_lijst <- resultaat_C57BL_2$log2FoldChange
names(gen_lijst) <- gsub("gene-", "",row.names(resultaat_C57BL_2))
gen_lijst <- sort(gen_lijst, decreasing = TRUE)
head(gen_lijst)## Myl3 Myh7 Fgb Fgg Csn3
## 24.57719 24.02822 23.08704 22.98055 22.10822
## 2310065F04Rik
## 20.36657Hier zie je een genlijst met genen met een log2FoldChange van hoog naar laag gesorteerd.
Met de genlijst doen we nu een gseGo met een cutofvalue van 0.05. Dat is normaal de p cut of waarde. En niet te vergeten dat we de muis gebruik Dus Mm en niet Hs gebruiken.
gsa_C57BL <- gseGO(geneList = gen_lijst,
keyType = "SYMBOL",
OrgDb = org.Mm.eg.db,
ont ="BP",
pvalueCutoff = 0.05,
verbose = T)
head(as.data.frame(gsa_C57BL))## ID Description setSize
## GO:0030216 GO:0030216 keratinocyte differentiation 152
## GO:0031424 GO:0031424 keratinization 59
## GO:0033561 GO:0033561 regulation of water loss via skin 40
## GO:0061436 GO:0061436 establishment of skin barrier 35
## GO:0045104 GO:0045104 intermediate filament cytoskeleton organization 90
## GO:0045109 GO:0045109 intermediate filament organization 71
## enrichmentScore NES pvalue p.adjust qvalue rank
## GO:0030216 -0.7074820 -3.520716 1e-10 1.899697e-08 1.485167e-08 2451
## GO:0031424 -0.8460711 -3.479936 1e-10 1.899697e-08 1.485167e-08 1145
## GO:0033561 -0.7897880 -2.968105 1e-10 1.899697e-08 1.485167e-08 1719
## GO:0061436 -0.7973078 -2.895061 1e-10 1.899697e-08 1.485167e-08 1485
## GO:0045104 -0.6278199 -2.820374 1e-10 1.899697e-08 1.485167e-08 2426
## GO:0045109 -0.6598070 -2.798259 1e-10 1.899697e-08 1.485167e-08 2426
## leading_edge
## GO:0030216 tags=47%, list=9%, signal=43%
## GO:0031424 tags=64%, list=4%, signal=62%
## GO:0033561 tags=52%, list=6%, signal=49%
## GO:0061436 tags=51%, list=5%, signal=49%
## GO:0045104 tags=43%, list=9%, signal=40%
## GO:0045109 tags=55%, list=9%, signal=50%
## core_enrichment
## GO:0030216 Sprr2d/Krt77/Csta2/Pou2f3/Notch1/Krt36/Krt8/Pou3f1/Vdr/Bcl11b/Foxc1/St14/Epha2/Cdh3/Irf6/Sfn/Trp63/Krt10/Tfap2c/Krt2/Ivl/Kdf1/Tfap2a/Dsg4/Tmem79/Krt7/Krt5/Sprr4/Tgm1/Stfa3/Grhl2/Sprr2f/Sprr2h/Krt90/Sprr2e/Scel/Krt14/Sprr2k/Krt6a/Krt80/Grhl1/Dsp/Stfa2/Gm5478/Krt79/Exph5/Foxn1/Cstdc5/Cers3/Csta1/Sprr1b/Evpl/Krt84/Gm5414/Ppl/Sprr2b/Sprr1a/Krt4/Krt1/Krt75/Cstdc6/Krt78/Cnfn/Krt16/Cdsn/Krt6b/Pnpla1/Sprr2i/Casp14/Sprr3/Tgm3/Abca12
## GO:0031424 Cdh3/Sfn/Krt2/Tmem79/Krt7/Krt5/Sprr4/Tgm1/Sprr2f/Sprr2h/Krt90/Sprr2e/Sprr2k/Krt6a/Krt80/Gm5478/Krt79/Cers3/Sprr1b/Evpl/Krt84/Gm5414/Ppl/Sprr2b/Sprr1a/Krt4/Krt1/Krt75/Krt78/Cnfn/Krt16/Cdsn/Krt6b/Sprr2i/Casp14/Sprr3/Tgm3/Abca12
## GO:0033561 Cdh1/Flg2/Alox12/Cldn4/Klf4/Sfn/Trp63/Kdf1/Tmem79/Tmprss13/Grhl1/Aloxe3/Grhl3/Krt1/Krt16/Acer1/Tmprss11f/Pnpla1/Flg/Alox12b/Abca12
## GO:0061436 Flg2/Alox12/Cldn4/Klf4/Sfn/Trp63/Kdf1/Tmem79/Tmprss13/Grhl1/Aloxe3/Grhl3/Krt1/Krt16/Pnpla1/Flg/Alox12b/Abca12
## GO:0045104 Krt77/Nefh/Krt33a/Krt36/Krt8/Tchh/BC024139/Krt19/Eppk1/Krt10/Krt2/Krt32/Krt40/Krt7/Krt5/Fam83h/Pkp1/Krt31/Krt90/Krt14/Krt6a/Krt80/Dsp/Gm5478/Krt15/Krt79/Evpl/Krt84/Gm5414/Ppl/Krt13/Krt4/Krt1/Krt75/Krt78/Krt23/Krt16/Krt6b/Flg
## GO:0045109 Krt87/Krt39/Krt33b/Krt74/Krt77/Nefh/Krt33a/Krt36/Krt8/Tchh/Krt19/Eppk1/Krt10/Krt2/Krt32/Krt40/Krt7/Krt5/Pkp1/Krt31/Krt90/Krt14/Krt6a/Krt80/Dsp/Gm5478/Krt15/Krt79/Krt84/Gm5414/Krt13/Krt4/Krt1/Krt75/Krt78/Krt23/Krt16/Krt6b/FlgHier zie je dus een Description met daaraan gekoppelde values die we gaan gebruiken.
Maar eerst even laten zien waar alles een beetje bij hoort.
boompie <- pairwise_termsim(gsa_C57BL)
treeplot(boompie)
Hier maken we een dotplot waar je kan zien van 5 Description of ze worden onderdrukt of juist meer tot uiting komen met de setSize en de generatio.
dotplot(gsa_C57BL, split=".sign", showCategory = 5) + facet_grid(.~.sign)
Hier zie je dus dat keratinization onderdrukt wordt. En de genen er voor veel voorkomen. En zie je ook dat striated muscle contraction meer tot uiting komt.
Hieronder laat ik zien welke genen er overeenkomen tussen keratinization en regulation of water loss via skin.
cnetplot(gsa_C57BL,
node_label="category",
showCategory = c("keratinization", "regulation of water loss via skin"))
cnetplot(gsa_C57BL,
node_label="gene",
showCategory = c("keratinization", "regulation of water loss via skin"))
Hieruit kan je concluderen dat die genen misschien minder tot uiting komen.
Hieronder laat ik dan de overeen komsten zien van myeloid leukocyte migration muscle contraction.
cnetplot(gsa_C57BL,
node_label="category",
showCategory = c("myeloid leukocyte migration", "muscle contraction"))
cnetplot(gsa_C57BL,
node_label="gene",
showCategory = c("myeloid leukocyte migration", "muscle contraction"))
Hieruit kan je concluderen dat die genen misschien meer tot uiting komen.
Hier maak ik een plot gsa van de bovenste 3 disciptions.
gseaplot2(gsa_C57BL, geneSetID = 1:3, subplots = 1:2)
In deze plot kan je zien waar de genen zich begeven in de genlijst en waar de es zich begeeft. De ES begeeft zich op het meest hoogste punt of net als in dit geval het laagste getal. Hier kan je dus zien dat keratinization de grootste afwijking heeft van de drie. De lijntjes eronder geven aan waar de genen zich begeven in de genlijst. Dus linksboven in de genlijst en rechtsonder in de genlijst die gesorteerd is op log2FoldChange.
Volcano
Hier zie je een eerste volcano plot van NSG_bladder tumor_Vehicle tegen C57BL/6_bladder tumor_Baseline. Dit had ik uit enthousiasme al gemaakt zonder te kijken wat er echt in de data stond.
EnhancedVolcano(res,
lab = rownames(res),
x = 'log2FoldChange',
y = 'pvalue',
title = "behandeld tegen onbehandeld",
pCutoff = 0.05,
FCcutoff = 1.5,
pointSize = 3.0,
labSize = 4.0,)
Hier zie je dus dat gene- Eno1b een grote invloed heeft, maar een mindere verandering heeft dan gen-Csn3.
Hieronder kijk ik welke hier de hoogste p-waarde hebben en zet ik die in de grafiek. Ik neem de 20 hoogste genen.
resultaat_oder <- res[order(res$pvalue), ]
top20genen <- rownames(head(resultaat_oder, 20))
EnhancedVolcano(res,
lab = rownames(res),
x = 'log2FoldChange',
y = 'pvalue',
selectLab = top20genen,
title = "behandeld tegen onbehandeld",
pCutoff = 0.05,
pointSize = 2.5,
labSize = 4
)
Hier kan je dus in 1 keer zien wat de stippen zijn met de hoogste p-waarde en kan je makkelijker een conclusie trekken.
Hieronder wordt de goede vergelijking gedaan. Dit zou ik moeten uitwerken.
EnhancedVolcano(resultaat_C57BL,
lab = rownames(resultaat_C57BL),
x = 'log2FoldChange',
y = 'pvalue',
title = "behandeld tegen onbehandeld C57BL ",
pCutoff = 0.05,
FCcutoff = 1.5,
pointSize = 2.5,
labSize = 4.0,
col=c('black', 'black', 'black', 'orange'))
Hier zie je dus dat er een grote verandering is bij het gen Myh7, dus is deze interessant om te bekijken. Myh7 gen is een gen codeert voor de belangrikjste hartspier. Dit geldt ook voor Bard1, want die heeft een hogere invloed erop. Bard1 gen heeft invoeld op gen respresie dus dat het gen meer tot uiting komt is een goed voordeel voor de het bestrijden van kanker.
MA
Hier ga ik en MA plot maken, zodat je dan kan zien wat de gemiddelde expressie niveaus zijn tegenover de expressie van de genen tussen de twee situaties.
p3 <-plotMA(resultaat_C57BL)
Deze plot ziet er nogal raar uit dus ik ga even kijken wat er anders kan.
Ik ben erachter gekomen dat je de data eerst moet Shrinken, zodat je consistent een betere uitput krijgt.
lfcschrink_subset <- lfcShrink(C57BL_subset, coef="treatment_Entinostat_vs_Baseline", type="apeglm")
p4 <-plotMA(lfcschrink_subset )
Dit ziet er beter uit en nu kan je concluderen dat er best wel een groot verschil is tussen de genen die tot uiting komen bij de verschillende omstandigheden.
Ik ga hier nog proberen om het mooier te laten zien met ggplot.
data_frame_resultaat_C57B <- as.data.frame(resultaat_C57BL)
data_frame_resultaat_C57B <- data_frame_resultaat_C57B[!is.na(data_frame_resultaat_C57B$padj),]
data_frame_resultaat_C57B$significant <- data_frame_resultaat_C57B$padj < 0.5
ggplot(data_frame_resultaat_C57B,
aes(x=baseMean, y=log2FoldChange))+
geom_point(aes(color = significant), alpha = 0.4, size = 1.5)+
scale_x_log10()+
geom_hline(yintercept = 0, color = "blue", linetype = "dashed")+
labs(title = 'maplot',
x = 'gemidelde',
y = "log2fold")+
scale_color_manual(values = c("TRUE" = "red", "FALSE" = "black"))+
theme_minimal()
Het verschil is er niet echt en het is naar mijn idee ook wel minder duidelijk. Want nu lijkt het alsof er best wel weinig niet anders is maar bij de MA plot functie juist wel te zien is.
Heatmap
Hier ga ik een heat map maken.
library(pheatmap)
topgene <- head(order(resultaat_C57BL$padj),20)
#mat <- assay(dds)[topgene,]
mat_2 <- assay(C57BL_subset)[topgene,]
#pheatmap(mat, cluster_rows=TRUE, cluster_cols = TRUE)
pheatmap(mat_2, cluster_rows=TRUE, cluster_cols = TRUE)
Deze vind ik nog al lelijk. Ik ga eerst even kijken naar een andere.
Ik heb een andere manier gevonden.
In de Heatpolt kan je ook makkelijk de Disciption meegeven die je graag zou willen zien, dus laat ik hier weer keratinization en regulation of water loss via skin zien.
heatplot(gsa_C57BL, foldChange=gen_lijst, showCategory=c("keratinization", "regulation of water loss via skin" ))
Hier kan je dus bij elk gen de invloed op de gegeven Disciption zien. In dit geval allemaal negatieve ### PCA
Hier maak ik een pca plot om te kijken of onze data een beetje goed verdeeld is.
vsd_C57BL <- vst(C57BL_subset,blind = FALSE)
DESeq2::plotPCA(vsd_C57BL, intgroup = "condition")
Bij de eerste oogopslag lijkt het goed, todat je je bedenkt dat er van de Baseline maar 4 zijn en van de entinostat 6. Hierdoor is het gemiddelde van de entinostat een stuk meer naar het midden dan van de Baseline.
Annotatie met biomart
Wij wouden gaan annoteren, omdat het ons interessant leek, maar na
overleg met de docenten hebben we besloten om het niet te doen. Ook
omdat het niet werkte, dus de code uit de meegegeven tutorial laat ik
wel staan, maar we doen er niks meer mee.
col_data <- col_data %>%
dplyr::group_by(strain, treatment) %>%
dplyr::mutate(r_num = row_number()) %>%
dplyr::ungroup() %>%
mutate(condition = paste0(strain, "_", treatment, "_r", r_num))
head(col_data)## # A tibble: 6 × 6
## Run source_name strain treatment r_num condition
## <chr> <chr> <chr> <chr> <int> <chr>
## 1 SRR12129014_star_ReadsPerGene C57BL/6_bladde… C57BL… Vehicle 1 C57BL/6_…
## 2 SRR12129015_star_ReadsPerGene C57BL/6_bladde… C57BL… Entinost… 1 C57BL/6_…
## 3 SRR12129016_star_ReadsPerGene C57BL/6_bladde… C57BL… Vehicle 2 C57BL/6_…
## 4 SRR12129017_star_ReadsPerGene C57BL/6_bladde… C57BL… Vehicle 3 C57BL/6_…
## 5 SRR12129018_star_ReadsPerGene C57BL/6_bladde… C57BL… Vehicle 4 C57BL/6_…
## 6 SRR12129019_star_ReadsPerGene C57BL/6_bladde… C57BL… Entinost… 2 C57BL/6_…for (i in 1:nrow(col_data)) {
idx <- grep(col_data$Run[i], names(counts))
names(counts)[idx] <- col_data$condition[i]
}
head(counts)## C57BL/6_Vehicle_r1 C57BL/6_Entinostat_r1 C57BL/6_Vehicle_r2
## MissingGeneID 5859 4387 7820
## gene-0610005C13Rik 3 0 0
## gene-0610006L08Rik 0 1 0
## gene-0610009E02Rik 0 1 1
## gene-0610009L18Rik 53 29 25
## gene-0610010K14Rik 1024 790 881
## C57BL/6_Vehicle_r3 C57BL/6_Vehicle_r4 C57BL/6_Entinostat_r2
## MissingGeneID 8707 3405 3821
## gene-0610005C13Rik 0 0 0
## gene-0610006L08Rik 0 0 2
## gene-0610009E02Rik 1 0 1
## gene-0610009L18Rik 93 43 34
## gene-0610010K14Rik 1408 897 523
## C57BL/6_Entinostat_r3 C57BL/6_Vehicle_r5 C57BL/6_Baseline_r1
## MissingGeneID 3777 6823 4673
## gene-0610005C13Rik 1 0 0
## gene-0610006L08Rik 0 0 0
## gene-0610009E02Rik 1 4 0
## gene-0610009L18Rik 44 93 25
## gene-0610010K14Rik 492 1304 975
## C57BL/6_Baseline_r2 C57BL/6_Entinostat_r4
## MissingGeneID 7162 2160
## gene-0610005C13Rik 0 1
## gene-0610006L08Rik 0 1
## gene-0610009E02Rik 0 1
## gene-0610009L18Rik 55 15
## gene-0610010K14Rik 1390 419
## C57BL/6_Baseline_r3 C57BL/6_Entinostat_r5
## MissingGeneID 5096 2718
## gene-0610005C13Rik 0 0
## gene-0610006L08Rik 0 0
## gene-0610009E02Rik 0 4
## gene-0610009L18Rik 24 49
## gene-0610010K14Rik 972 676
## C57BL/6_Entinostat_r6 C57BL/6_Vehicle_r6 C57BL/6_Baseline_r4
## MissingGeneID 3870 9562 5163
## gene-0610005C13Rik 0 0 0
## gene-0610006L08Rik 0 0 0
## gene-0610009E02Rik 6 0 2
## gene-0610009L18Rik 34 43 45
## gene-0610010K14Rik 431 964 876
## NSG_Baseline_r1 NSG_Baseline_r2 NSG_Baseline_r3
## MissingGeneID 8453 8045 5803
## gene-0610005C13Rik 4 1 2
## gene-0610006L08Rik 0 0 0
## gene-0610009E02Rik 7 4 7
## gene-0610009L18Rik 38 47 31
## gene-0610010K14Rik 967 1298 664
## NSG_Baseline_r4 NSG_Baseline_r5 NSG_Entinostat_r1
## MissingGeneID 6777 8639 6767
## gene-0610005C13Rik 1 3 1
## gene-0610006L08Rik 1 0 0
## gene-0610009E02Rik 8 10 6
## gene-0610009L18Rik 45 48 95
## gene-0610010K14Rik 1001 1336 1470
## NSG_Entinostat_r2 NSG_Entinostat_r3 NSG_Entinostat_r4
## MissingGeneID 9366 9204 3295
## gene-0610005C13Rik 6 2 0
## gene-0610006L08Rik 0 0 0
## gene-0610009E02Rik 5 6 0
## gene-0610009L18Rik 139 124 64
## gene-0610010K14Rik 1528 1339 387
## NSG_Entinostat_r5 NSG_Vehicle_r1 NSG_Vehicle_r2
## MissingGeneID 12470 5393 7122
## gene-0610005C13Rik 2 2 1
## gene-0610006L08Rik 1 0 0
## gene-0610009E02Rik 13 6 11
## gene-0610009L18Rik 113 41 50
## gene-0610010K14Rik 1393 735 1145
## NSG_Vehicle_r3 NSG_Vehicle_r4 NSG_Vehicle_r5
## MissingGeneID 8495 8854 7730
## gene-0610005C13Rik 2 4 2
## gene-0610006L08Rik 0 0 1
## gene-0610009E02Rik 5 8 4
## gene-0610009L18Rik 122 150 160
## gene-0610010K14Rik 1618 1689 1511counts$Ensembl <- mapIds(x = org.Mm.eg.db,
keys=gsub("gene-", "",row.names(counts)),
column="ENSEMBL",
keytype="SYMBOL",
multiVals="first")library(biomaRt)
ensembl=useMart("ENSEMBL_MART_ENSEMBL", host="https://www.ensembl.org")
ensembl <- useMart("ensembl")
mart.datasets <- listDatasets(ensembl)
ensembl <- useDataset('mmusculus_gene_ensembl', mart = ensembl)
filters <- listFilters(ensembl)
attributes <- listAttributes(ensembl)# Set the 'attributes' values
attrs.get <- c("ensembl_gene_id", "chromosome_name",
"start_position","end_position", "description")
Perform a biomaRt query using 'getBM'
results <- getBM(attributes = attrs.get,
filters = "ensembl_gene_id",
values = counts$Ensembl[12:15],
mart = ensembl, verbose = TRUE)
results$gene_length <- abs(results$end_position - results$start_position)
merge(x = counts, y = results, by.x = 'Ensembl', by.y = 'ensembl_gene_id')